home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / listings / v_13_11 / phillip2 / fitt.c < prev    next >
C/C++ Source or Header  |  1993-08-20  |  18KB  |  650 lines
  1. /* (C) Copr. 1986-92 Numerical Recipes Software  */
  2.  
  3.  
  4.    /********************************************
  5.    *
  6.    *   This file contains code from the 
  7.    *   Numerical Recipies in C text by
  8.    *   Press, Flannery, Teukolsky, Vetterling
  9.    *   I needed the routines to fit data to
  10.    *   a straight line y=mx+b.
  11.    *
  12.    *
  13.    *   This files contains the files:
  14.    *      nrutil.h
  15.    *      fit.c
  16.    *      gammq.c
  17.    *      gcf.c
  18.    *      gser.c
  19.    *      gammln.c
  20.    *      nrutil.c
  21.    *
  22.    *   21 August 1993
  23.    *
  24.    ********************************************/
  25.  
  26.  
  27.  
  28.  
  29. #include <math.h>
  30.  
  31.    /* file nrutil.h */
  32.  
  33. #ifndef _NR_UTILS_H_
  34. #define _NR_UTILS_H_
  35.  
  36. static float sqrarg;
  37. #define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)
  38.  
  39. static double dsqrarg;
  40. #define DSQR(a) ((dsqrarg=(a)) == 0.0 ? 0.0 : dsqrarg*dsqrarg)
  41.  
  42. static double dmaxarg1,dmaxarg2;
  43. #define DMAX(a,b) (dmaxarg1=(a),dmaxarg2=(b),(dmaxarg1) > (dmaxarg2) ?\
  44.         (dmaxarg1) : (dmaxarg2))
  45.  
  46. static double dminarg1,dminarg2;
  47. #define DMIN(a,b) (dminarg1=(a),dminarg2=(b),(dminarg1) < (dminarg2) ?\
  48.         (dminarg1) : (dminarg2))
  49.  
  50. static float maxarg1,maxarg2;
  51. #define FMAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > (maxarg2) ?\
  52.         (maxarg1) : (maxarg2))
  53.  
  54. static float minarg1,minarg2;
  55. #define FMIN(a,b) (minarg1=(a),minarg2=(b),(minarg1) < (minarg2) ?\
  56.         (minarg1) : (minarg2))
  57.  
  58. static long lmaxarg1,lmaxarg2;
  59. #define LMAX(a,b) (lmaxarg1=(a),lmaxarg2=(b),(lmaxarg1) > (lmaxarg2) ?\
  60.         (lmaxarg1) : (lmaxarg2))
  61.  
  62. static long lminarg1,lminarg2;
  63. #define LMIN(a,b) (lminarg1=(a),lminarg2=(b),(lminarg1) < (lminarg2) ?\
  64.         (lminarg1) : (lminarg2))
  65.  
  66. static int imaxarg1,imaxarg2;
  67. #define IMAX(a,b) (imaxarg1=(a),imaxarg2=(b),(imaxarg1) > (imaxarg2) ?\
  68.         (imaxarg1) : (imaxarg2))
  69.  
  70. static int iminarg1,iminarg2;
  71. #define IMIN(a,b) (iminarg1=(a),iminarg2=(b),(iminarg1) < (iminarg2) ?\
  72.         (iminarg1) : (iminarg2))
  73.  
  74. #define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
  75.  
  76. #if defined(__STDC__) || defined(ANSI) || defined(NRANSI) /* ANSI */
  77.  
  78. void nrerror(char error_text[]);
  79. float *vector(long nl, long nh);
  80. int *ivector(long nl, long nh);
  81. unsigned char *cvector(long nl, long nh);
  82. unsigned long *lvector(long nl, long nh);
  83. double *dvector(long nl, long nh);
  84. float **matrix(long nrl, long nrh, long ncl, long nch);
  85. double **dmatrix(long nrl, long nrh, long ncl, long nch);
  86. int **imatrix(long nrl, long nrh, long ncl, long nch);
  87. float **submatrix(float **a, long oldrl, long oldrh, long oldcl, long oldch,
  88.    long newrl, long newcl);
  89. float **convert_matrix(float *a, long nrl, long nrh, long ncl, long nch);
  90. float ***f3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh);
  91. void free_vector(float *v, long nl, long nh);
  92. void free_ivector(int *v, long nl, long nh);
  93. void free_cvector(unsigned char *v, long nl, long nh);
  94. void free_lvector(unsigned long *v, long nl, long nh);
  95. void free_dvector(double *v, long nl, long nh);
  96. void free_matrix(float **m, long nrl, long nrh, long ncl, long nch);
  97. void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch);
  98. void free_imatrix(int **m, long nrl, long nrh, long ncl, long nch);
  99. void free_submatrix(float **b, long nrl, long nrh, long ncl, long nch);
  100. void free_convert_matrix(float **b, long nrl, long nrh, long ncl, long nch);
  101. void free_f3tensor(float ***t, long nrl, long nrh, long ncl, long nch,
  102.    long ndl, long ndh);
  103.  
  104. #else /* ANSI */
  105. /* traditional - K&R */
  106.  
  107. void nrerror();
  108. float *vector();
  109. float **matrix();
  110. float **submatrix();
  111. float **convert_matrix();
  112. float ***f3tensor();
  113. double *dvector();
  114. double **dmatrix();
  115. int *ivector();
  116. int **imatrix();
  117. unsigned char *cvector();
  118. unsigned long *lvector();
  119. void free_vector();
  120. void free_dvector();
  121. void free_ivector();
  122. void free_cvector();
  123. void free_lvector();
  124. void free_matrix();
  125. void free_submatrix();
  126. void free_convert_matrix();
  127. void free_dmatrix();
  128. void free_imatrix();
  129. void free_f3tensor();
  130.  
  131. #endif /* ANSI */
  132.  
  133. #endif /* _NR_UTILS_H_ */
  134.  
  135.  
  136.  
  137.  
  138.  
  139.    /* file fit.c */
  140.  
  141.  
  142. #undef  DDEBUG
  143. #define NRANSI
  144.  
  145. void fit(float x[], float y[], int ndata, float sig[], int mwt, float *a,
  146.    float *b, float *siga, float *sigb, float *chi2, float *q)
  147. {
  148.    float gammq(float a, float x);
  149.    int i;
  150.    float wt,t,sxoss,sx=0.0,sy=0.0,st2=0.0,ss,sigdat;
  151.  
  152. int iii;
  153. #ifdef DDEBUG
  154. printf("\nDEBUG> ndata=%d mwt=%d", ndata, mwt);
  155. for(iii=1; iii<=ndata; iii++)
  156. printf("\nDEBUG> x=%f y=%f sig=%f",x[iii],y[iii],sig[iii]);
  157. printf("\nDEBUG> start of fit");
  158. #endif DDEBUG
  159.  
  160. /* 
  161.    change 8-21-93
  162.    The sig array is all 1's so I will remove
  163.    all uses of it in mult and division.
  164. */
  165.  
  166.    *b=0.0;
  167.    if (mwt) {
  168.       ss=0.0;
  169.       for (i=1;i<=ndata;i++) {
  170.          wt=1.0;
  171.          ss += wt;
  172.          sx += x[i]*wt;
  173.          sy += y[i]*wt;
  174.       }
  175.    } else {
  176.       for (i=1;i<=ndata;i++) {
  177.          sx += x[i];
  178.          sy += y[i];
  179.       }
  180.       ss=ndata;
  181.    }
  182.    sxoss=sx/ss;
  183.    if (mwt) {
  184.       for (i=1;i<=ndata;i++) {
  185.          t=(x[i]-sxoss);
  186.          st2 += t*t;
  187.          *b += t*y[i];
  188.       }
  189.    } else {
  190.       for (i=1;i<=ndata;i++) {
  191.          t=x[i]-sxoss;
  192.          st2 += t*t;
  193.          *b += t*y[i];
  194.       }
  195.    }
  196.    *b /= st2;
  197.    *a=(sy-sx*(*b))/ss;
  198.    *siga=sqrt((1.0+sx*sx/(ss*st2))/ss);
  199.    *sigb=sqrt(1.0/st2);
  200.    *chi2=0.0;
  201.    if (mwt == 0) {
  202.       for (i=1;i<=ndata;i++)
  203.          *chi2 += SQR(y[i]-(*a)-(*b)*x[i]);
  204.       *q=1.0;
  205.       sigdat=sqrt((*chi2)/(ndata-2));
  206.       *siga *= sigdat;
  207.       *sigb *= sigdat;
  208.    } else {
  209.       for (i=1;i<=ndata;i++)
  210.          *chi2 += SQR((y[i]-(*a)-(*b)*x[i]));
  211.       *q=gammq(0.5*(ndata-2),0.5*(*chi2));
  212.    }
  213. }
  214. #undef NRANSI
  215.  
  216.  
  217.  
  218.  
  219.  
  220.    /* file gammq.c */
  221.  
  222. /* (C) Copr. 1986-92 Numerical Recipes Software $|9`0S[)03. */
  223. float gammq(float a, float x)
  224. {
  225.    void gcf(float *gammcf, float a, float x, float *gln);
  226.    void gser(float *gamser, float a, float x, float *gln);
  227.    void nrerror(char error_text[]);
  228.    float gamser,gammcf,gln;
  229.  
  230.    if (x < 0.0 || a <= 0.0) nrerror("Invalid arguments in routine gammq");
  231.    if (x < (a+1.0)) {
  232.       gser(&gamser,a,x,&gln);
  233.       return 1.0-gamser;
  234.    } else {
  235.       gcf(&gammcf,a,x,&gln);
  236.       return gammcf;
  237.    }
  238. }
  239.  
  240.  
  241.  
  242.  
  243.  
  244.    /* file gcf.c */
  245.  
  246. /* (C) Copr. 1986-92 Numerical Recipes Software $|9`0S[)03. */
  247. #define ITMAX 100
  248. #define EPS 3.0e-7
  249. #define FPMIN 1.0e-30
  250.  
  251. void gcf(float *gammcf, float a, float x, float *gln)
  252. {
  253.    float gammln(float xx);
  254.    void nrerror(char error_text[]);
  255.    int i;
  256.    float an,b,c,d,del,h;
  257.  
  258.    *gln=gammln(a);
  259.    b=x+1.0-a;
  260.    c=1.0/FPMIN;
  261.    d=1.0/b;
  262.    h=d;
  263.    for (i=1;i<=ITMAX;i++) {
  264.       an = -i*(i-a);
  265.       b += 2.0;
  266.       d=an*d+b;
  267.       if (fabs(d) < FPMIN) d=FPMIN;
  268.       c=b+an/c;
  269.       if (fabs(c) < FPMIN) c=FPMIN;
  270.       d=1.0/d;
  271.       del=d*c;
  272.       h *= del;
  273.       if (fabs(del-1.0) < EPS) break;
  274.    }
  275.    if (i > ITMAX) nrerror("a too large, ITMAX too small in gcf");
  276.    *gammcf=exp(-x+a*log(x)-(*gln))*h;
  277. }
  278. #undef ITMAX
  279. #undef EPS
  280. #undef FPMIN
  281.  
  282.  
  283.  
  284.  
  285.  
  286.    /* file gser.c */
  287.  
  288. /* (C) Copr. 1986-92 Numerical Recipes Software $|9`0S[)03. */
  289. #define ITMAX 100
  290. #define EPS 3.0e-7
  291.  
  292. void gser(float *gamser, float a, float x, float *gln)
  293. {
  294.    float gammln(float xx);
  295.    void nrerror(char error_text[]);
  296.    int n;
  297.    float sum,del,ap;
  298.  
  299.    *gln=gammln(a);
  300.    if (x <= 0.0) {
  301.       if (x < 0.0) nrerror("x less than 0 in routine gser");
  302.       *gamser=0.0;
  303.       return;
  304.    } else {
  305.       ap=a;
  306.       del=sum=1.0/a;
  307.       for (n=1;n<=ITMAX;n++) {
  308.          ++ap;
  309.          del *= x/ap;
  310.          sum += del;
  311.          if (fabs(del) < fabs(sum)*EPS) {
  312.             *gamser=sum*exp(-x+a*log(x)-(*gln));
  313.             return;
  314.          }
  315.       }
  316.       nrerror("a too large, ITMAX too small in routine gser");
  317.       return;
  318.    }
  319. }
  320. #undef ITMAX
  321. #undef EPS
  322.  
  323.  
  324.  
  325.  
  326.  
  327.    /* file gammln.c */
  328.  
  329. /* (C) Copr. 1986-92 Numerical Recipes Software $|9`0S[)03. */
  330.  
  331. float gammln(float xx)
  332. {
  333.    double x,y,tmp,ser;
  334.    static double cof[6]={76.18009172947146,-86.50532032941677,
  335.       24.01409824083